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I . Introduction 

This report contains user instructions, a listing and 
documentation for a microcomputer BASIC program that can be used 
to compute an estimate of the probability that a magnetic 
anomaly detection (MAD) system such as the AN/ASQ-81 will detect 
a submarine during an encounter. 

The program generates detection probabilities based on two 
encounter models. In the first encounter model, the detection 
system uses a crosscorrelation detector. In the second 
encounter model, the detection system uses a square law detector. 
Relative to operationally realizable values, probabilities based 
on the first model represent upper bounds and those based on the 
second model represent lower bounds. For both encounter models, 
a magnetometer signal is proportional to the magnitude of the 
projection of an anomaly field on the earth field, a submarine 
anomaly field is a dipole field and the earth field and the noise 
do not change with changes in the position of the magnetometer. 
Also, for both encounter models, an encounter is a straight line 
encounter with constant vertical separation. 

The encounter models can be interpreted as models of a 
magnetic anomaly detection system on an aircraft that is moving 
with constant course, speed and altitude in an encounter with a 
submarine moving with constant course, speed and depth. Or, they 
can be interpreted as models of a stationary magnetic anomaly 
detection system in an encounter with a submarine moving with 
constant course, speed and depth. 
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The program parameters include encounter latitude and 
longitude, submarine induced magnetic moments, submarine 
permanent magnetic moments, submarine course, speed and altitude, 
magnetometer course, speed and depth, encounter lateral range 
(the horizontal range at the closest point of approach in a 
straight line encounter) and false alarm rate. (In the program, 
a false alarm is the event that the detection system classifies 
noise as a dipole signal.) 
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II. Program User Instructions 

As listed in Appendix 10, the program can be run under a 
BASIC language that is compatible with IBM PC BASIC. If the 
listing is used to enter the program through the keyboard, then 
the program should be saved with the name MAD. BAS or the value of 
N$ on line 20 should be changed to the file name under which it 
is saved. 

The program contains user instructions in the form of query 
and parameter limitation messages. As an example of the former, 
after starting the program under BASIC, the following message 
should appear: 

Magnetic Anomaly Detection (MAD) Lateral Range Function 

generate or print a program data file (g/p)? 

By entering p, data can be printed from a program data file 
that was generated by the program. By entering g, a program 
data file can be generated for a set of user specified 
conditions. With either response, a sequence of additional 
queries is displayed. These queries require either an indicated 
response or a parameter value as the input. If the initial 
response is g, the sequence includes queries whose responses 
determine whether or not an auxiliary data file will be generated 
that can be used for future input of magnetic, processing or 
kinematic data. In particular, the first query in the sequence 
gives the option of using a combined magnetic, processing and 
kinematic data file. If the response indicates that it should 
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be used and the file is available, the parameter values that 
remain to be entered in order to generate a program data file are 
the following: the false alarm rate, the magnetic noise, the 
maximum encounter lateral range and the lateral range step. The 
combined file should be used only if the effect of varying just 
one or more of these parameters is desired. If the response 
indicates that the file will not be used, queries concerning 
magnetic, kinematic and processing parameter values are 
displayed. 

After all of the program parameter values have been 
entered, a query is displayed giving the option of generating the 
combined file. Then a query giving the option of generating a 
program data file, a query giving the option of printing the 
encounter parameter values and a query giving the option of 
printing lateral range function values are displayed. The 
lateral range function values are the encounter detection 
probabilities indexed by lateral range. The parameters maximum 
lateral range and lateral range step determine the index lateral 
ranges of the encounters for which probabilities of detection are 
computed . 

The program generates magnetic signal values that correspond 
to points in time during an encounter. Following the lateral 
range function query, a query is displayed that gives the option 
of printing the magnetic signal values for an encounter. If the 
option is exercised, the option is repeated. When the option is 
not exercised, a query is displayed that gives the option of 
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generating or printing a new program data file. If this option 
is not exercised, the program ends. 

Some suggested guides for determining parameter values can 
be found in Section IV of this report. 
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III. Encounter Model Limitations 



In the two encounter models, noise is determined by a 
gaussian process. This may be a limitation in describing a 
stationary magnetometer. In this case, the noise must account 
for geomagnetic, instrument and magnetohydrodynamic noise. And 
it may be a significant limitation if the magnetometer is 
moving. In this case, in addition to the above sources, the 
noise must account for maneuver noise and geologic noise. 

The signal is that of a magnetic dipole that moves relative 
to the magnetometer. Describing a ship magnetic anomaly as a 
dipole anomaly should not be a significant limitation for 
encounter slant range at the closest point of approach (CPA) that 
is greater than one hull length. The detection decision is 
based on samples from a single time interval (window) that is 
centered on the CPA. The length of the sampling interval and the 
sampling rate are parameter values that are inputs to the 
program. In terms of signal-to-noise ratio, there is an optimum 
sampling interval length (integration time) and sample rate. 
Although dipole signal energy is not symmetrically distributed 
about the CPA time, for a given sampling interval, the difference 
between the signal energy for the optimum interval location and 
the CPA centered location should not be significant in most 
cases . 

The encounter model magnetometer noise samples are values of 
independent identically distributed random variables. The 
standard deviation of these random variables is referred to as 
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the magnetic noise and the variance is the noise in the sense of 
the signal-to-noise ratio. The magnetic noise of the encounter 
models is determined by a gaussian random process. The degree of 
correspondence between this process and operational noise 
depends on the nature of the dominant operational noise sources 
and on the magnetometer input filtering (noise whitening) . 

In the encounter models, the intervals are adjacent but not 
overlapping and a detection statistic corresponds to each sample 
interval and its value is determined by the sample values. If 
the value of the detection statistic for an interval equals or 
exceeds a threshold value, a detection is indicated. The 
threshold value is determined by the false alarm probability 
which in turn is determined by the false alarm rate and the 
sampling interval length. 

A detector that used a moving sample interval ' that was 
generated by replacing the oldest sample by the newest one would 
correspond more closely to the detector in an operational 
detection system. In a model of the detector, since the sample 
windows overlap, the detection statistics would represent a 
sequence of random variables that were correlated over an 
interval equal to the width of the sample interval. Because of 
this dependence, it seems unlikely that the results that would be 
obtained with an encounter model based on an overlapping interval 
detector would differ significantly from those obtained with an 
encounter model based on a nonoverlapping interval detector. 
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For encounter lengths of the order of a few nautical miles 
or less, the straight line encounter condition should not be a 
significant limitation. In particular, this should be the case 
for a fixed magnetometer since, for a submarine (or surface ship) 
target, vertical separation and course changes should be less 
likely to occur. 

Other models are available that can be used as the basis for 
computing an estimate of the probability that a magnetic anomaly 
detection system will detect a submarine during an encounter. 

For example, one that is described in Appendix 8 can be used to 
determine the slant ranges of straight line encounters for which 
the detection probability is equal to a specified value. The 
parameter values that are required to do this are an average 
submarine dipole moment, a detection system capability factor and 
a noise factor. Values for these parameters can be determined 
from operational data. However, the values are specific to 
averages over a particular set of encounter conditions. An 
advantage of the two encounter models relative to this model is 
their adaptability to different magnetic, processing and 
kinematic conditions. 
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IV. Parameter Values 



The magnetic parameter value queries are generally 
explanatory with regard to the value that should be entered. 

This is also true of the kinematic parameter value queries. 
However, there is some ambiguity with respect the processing 
parameter value queries and the noise parameter value query. To 
reduce this ambiguity, a brief discussion of the common 
characteristics of the two encounter model processing parameters 
and the noise parameter is given below. This is followed by 
some guidelines for choosing these parameter values. 

In both encounter models, a decision is made at the end of 
each sampling interval. The decision is either noise energy was 
present during the interval or noise energy plus signal energy 
was present during the interval. The sample intervals are 
adjacent, equal width, nonoverlapping time intervals. The number 
of samples that are input in a sample interval is determined by 
the sampling rate and the interval length. 

The program default choice for the sampling rate is 2 ■ MAXF 
where MAXF is a parameter that is labeled the maximum magnetic 
signal frequency. This sample rate is the Nyquist rate for an 
ideal low pass filter. However, the signal in a sample interval 
that is computed by the program represents an unfiltered dipole 
signal. This is a reasonable approximation if the signal energy 
that is associated with signal components greater than MAXF is 
relatively small. As discussed below, the noise energy in a 
sample interval should be considered to be proportional to MAXF 
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in order to be consistent with the encounter models. Ideally, a 
default choice for MAXF would make the ratio of the signal energy 
to the noise energy a maximum for the sampling interval of an 
encounter. The program default choice for MAXF is 2 • MAXVM/MINRO 
where MAXVM is a user estimated maximum encounter relative speed 
converted from knots to meters per second and MINRO is a user 
estimated maximum slant range at CPA in meters in terms of a just 
detectable target. (A precise definition of a just detectable 
target can be made in terms of a specified detection probability, 
false alarm probability and target dipole moment.) The default 
choice for MAXF is consistent with the observation in Reference 1 
that if an optimum value for MAXF is determined for a minimum 
dipole moment target, then no significant increase in MAXF is 
required in order to maintain a required detection probability if 
the encounter lateral range is decreased even though the signal 
energy spectrum is shifted to higher frequencies. 

Because of the detection statistics and the gaussian noise 
model, if there were no penalty for decision delay, a sample 
interval length for a signal should be chosen equal to the signal 
duration, since this would make the detection probability for an 
encounter a maximum. In the program, a default choice for the 
sample interval length (integration time) is 2 -MAXRO/MINVM where 
MINVM is a user estimated maximum encounter relative speed 
converted from knots to meters per second and MAXRO is a user 
estimated maximum slant range at CPA in meters for a detectable 
encounter (in terms of a specific detection probability and false 
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alarm probability) with a user specified maximum dipole moment 
target. This choice might be considered a balance between 
minimizing decision delay and maximizing detection probability. 

In the encounter models, the sample intervals are located so that 
the CPA time is at the center of a sample interval and this 
sample interval is the only one that contains signal energy. 

This characteristic is consistent with the default choice for the 
sample interval. For a given sample interval length, although in 
general the interval location is not the optimum one in terms of 
signal energy, it should be approximately so in most cases. 

The program noise parameter is SIG. It represents the 
standard deviation a associated with the magnetic noise process 
of the two encounter models. In terms of the ideal low pass 
filter implied by the encounter models, its square should be 
equal to MAXF- (SIGO) * where SIGO is the magnetic noise process 
constant power spectral density. The program does not enforce 
this relation. Therefore, in using the program, the implied 
relation between the two input parameters: magnetic noise and 
maximum magnetic signal frequency should be kept in mind. If an 
average value of the peak-to-peak magnetic noise for an encounter 
can be estimated, for example from a magnetometer trace, then the 
value for SIG should be chosen so that the estimate is 4 to 6 
times this value. 
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Appendix 1. The Detection Statistics 

In this appendix, Yii Y2> ' ' ' • sequential values 

(voltages) representing the sample values in a sample interval. 
They are the input to a magnetometer's detector. With these m 
sample values, the detector computes the value of a detection 
statistic. This value is represented by x and the detection 
decision corresponding to the sample interval is determined by 
the decision rule: If x > x*, then the input during the sample 
interval was noise plus signal, otherwise, the input was noise. 

For both encounter models, the detection probability and the 
false alarm probability are decreasing functions of x* and the 
relation is one-to-one in both cases. In the program, the false 
alarm probability pf is used to determine a unique value of 
the threshold x*. This value is then used to determine a unique 
value of the detection probability p^j. 

In the program, Pf is found using the relation Pf = R-5t 
where R is the false alarm rate in false alarms per second 
and St is the sample interval length in seconds. This relation 
is based on the following argument: With no signal energy in a 
sample interval, yi = n^, Y2 = ri2 ' ' ‘ ‘ ' ^m “ ^m where n^^, n 2 , 

• • • , nj^ are noise values (voltages) input to a magnetometer's 

detector. The nj^ are values of independent normal (gaussian) 
random variables, each with mean zero and standard deviation a. 
Because of this, in the encounter models, values of x for 
different sample intervals are the values of independent random 
variables that determine two outcomes: x > x* or x < x*. 
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Therefore, in terms of these outcomes, an encounter is a sequence 
of independent Bernoulli trials. If there is no signal present, 
since the noise is determined by a stationary process, pf is the 
same for each sample interval and, in terms of these outcomes, 
the sequence is a series of repeated independent Bernoulli trials 
and therefore the expected number of trials between false alarms 
is 1/Pf* Since the time between trials is 5t, the expected 
number of seconds between false alarm is 5t/pf or the expected 
number of false alarms per second R is equal to pf/5t. 

The determination of x* depends on the encounter model 
statistic. For both encounter models, when there is a signal, 

yi = ni + si, y 2 = n 2 + y 2 / • • • > Ym = where s^, S 2 , 

••• , Sj(, are signal values (voltages) input to a magnetometer's 

detector. The models imply that the signal values Sj^ = K- (Hg) 
where K is a constant whose value is determined by the 
characteristics of the encounter magnetometer and where the 
(Hs) i are dipole magnetic signal intensities. The models also 
imply that the noise values nj_ are determined by a gaussian 
stochastic processes characterized by a standard deviation a 
and that nj_ = + n]_ where the (H^) i are magnetic noise 

intensities and the n| are magnetometer instrument noise values 
(voltages) . In the program, the magnitude of K is 1. Since, 
for both models, p^j depends only on the ratio of signal energy 
to noise energy for a sample interval, this is a satisfactory 
choice for the program if instrument noise is assumed to be 
determined by a process independent of the magnetic noise process 
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and to be expressed in terms of an equivalent magnetic noise 
(H ) [ = (1/K) -ni. 

For both encounter models, the signal (the average signal 

2 

power) S = (l/m)-S sj_ and the noise (the expected value of the 

2 

average noise power) N = a so that the signal-to-noise ratio 
IS (1/m) -S sj[/a where the sum index i = 1, 2, • • • , m. 

The Crosscorrelation Detector Statistic: The statistic for the 

first encounter model is a crosscorrelation detector statistic 
that is defined by the sum 

X = 2 Yi'Si 

where the summation index i = 1, 2, • • • , m and the sum is over 

the values corresponding to a sample interval. For the first 
encounter model , the characteristics of both the noise and the 
signal are required in order to determine encounter detection 
probabilities. In particular, the signal values for an encounter 
are in the memory of the detector prior to the encounter. For 
the encounter conditions and a specified false alarm 
probability, the statistic is optimum in the sense that the 
encounter detection probability for this statistic is at least 
equal to that for any other statistic. Because of these 
considerations, encounter probabilities based on the 
crosscorrelation statistic can be considered to represent upper 
bounds on detection performance against dipole targets for 
magnetometers of the type described by the models. 

For a sample interval without signal energy, x is the value 
of a normal random variable with a mean = 0 and a variance 
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2 2 2 

0 -^ = 0 -S Si where a is the standard deviation associated with 
the noise process and the sum index 1=1, 2, • • • ,m and the 

sum is over the values corresponding to the sample interval. 

This implies that 

Pf = 1 - P(x*/Ox) 

where P(z) is the standard normal cumulative distribution 
function. This relation is the basis for determining the 
threshold value x* . 

For the sample interval with signal energy, x is the value 

, , 2 

of a normal random variable with a mean ii^ = 2 s^ where the sum 
index i = 1, 2, • • • , m and the sum is over the values 

corresponding to the sample interval. This implies that 

Pd = 1 - P(v* - d^) 

where v* = x*/aj^ and d = S s^/a* = (1/m) -S/N. This relation 
is the basis for determining encounter detection probabilities 
for the first encounter model. The relation implies that for a 
specified false alarm probability pf the detection probability 
P(j is an increasing function of the signal to-noise ratio S/N. 
The Square Law Detector Statistic: The statistic for the second 

encounter model is a square law (energy) detector statistic that 
is defined by the sum 

X = S yi 

where the sum index i = 1, 2, • • • , m and the sum is over the 

values corresponding to the sample interval. For the second 
encounter model, only the characteristics of the noise are 
required to determine encounter detection probabilities. 
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For a sample interval without signal energy, x/a* is the 
value of a chi-square random variable with m degrees of 
freedom. This implies that 

Pf = 1 - P(x*/o" Im) 

where P(x /a |m) is the chi-square cumulative distribution 
function for a chi-square random variable with m degrees of 
freedom and where a is the standard deviation associated with 
the noise process. This relation is the basis for determining 

the threshold value x* . 

. . . 2 
For the sample interval with signal energy, x/a is the 

value of a noncentral chi-square random variable with m 

2 2 

degrees of freedom and noncentral parameter S s^/o where the 
sum index i = 1, 2, • • • , m and the sum is over the values 

corresponding to the sample interval. This implies that 

Pd = 1 - P(x*/cJ*|in^ 2 s\/a^) 

where P(x /a |m, S s^/o ) is the noncentral chi-square 

cumulative distribution function for a noncentral chi-square 

random variable with m degrees of freedom and noncentral 
2 2 

parameter S Sj_/a = (l/m)-S/N. This relation is the basis for 

determining encounter detection probabilities for the second 

encounter model. The relation implies that for a specified false 

alarm probability pf, the detection probability p^j is an 

increasing function of the signal-to-noise ratio S/N. This is 

made more evident by the following relation: 

P(x*/a"|m, 2 Si/a") = 2 { (aVj I ) ' exp (-a) • P[x*/a* 1 (m + 2-j)]} 

2 2 

where a = (1/2) -2 s^/a = (m/2) (S/N), the sum index i = 1, 2, 
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• • • , m and the sum index j = 0, 1, 2, ••• . (Note, > pj, 

as expected, since P(x*/a’|m) < P[x*/a’|(m + 2-j)] for j = 0, 
1 , 2 , • • • . ) 
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Appendix 2. Program Probability Calculations 

The program evaluates the cumulative and inverse cumulative 
distribution functions using approximations described in 
Reference 2. These approximations are listed below. 

The Standard Normal Cumulative Distribution Function 
Approximation ; 

P(z) = 1 - s-t-(bi + t- (b 2 + t- (b 3 + t- (b 4 + t-b 5 )))) 
where s = (1/(2 -tt) -exp(-z*/2) and t = 1/(1 + bQ-z). And 
where 

bo = .2316419, b]_ = .319381530, b2 = -.356563782, 

b4 = -1.821255978, bs = 1.330274429, 
and z > 0. For z < 0, P(z) = 1 - P(lz|). 

The Inverse Standard Normal Ciomulative Distribution Function 
Approximation : 

z(P) = t - (Co + t-(Ci + t-C2))/(l + t-(di + t-(d2 + td3))) 
where t = (ln(l/Q*))^ and Q = 1 - P. And where 
Co = 2.515517, Cl = .802853, C 2 = .010328, 
di = 1.432788, d2 = .189269, do = .001308, 
and .5 < P < 1. For 0 < P < .5, z(P) = -z(l - P) . 

The Inverse Chi-Square Ciimulative Distribution Function 
Approximation : 

v(p|m) = m-[l - 2/(9-m) + z • (2/ (9 • m) ) I ] ^ 
where P(z) = P(v|m). In the program, the inverse standard 
normal cumulative distribution function approximation is used to 
determine z . 
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The Noncentral Chi-Square Cumulative Distribution Function 
Approximation : 

P(w|m,S s\/o^) = P(z) 

where z = [2-w/(l + b)]i - [2-a/(l + b) - l]i with 

2 2 2 2 2 2 
a = m + S Sj_/a , b = (S s^/a )/(m + S s^/a ) and the sum index 

i=l, 2, In the program, the standard normal 

cumulative distribution function approximation is used to 

determine P(z). 
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Appendix 3 . The Magnetic Signal 

The encounter models are defined by the following 
conditions: A submarine magnetic anomaly field is a magnetic 

dipole field that is superimposed on an earth magnetic field that 
is constant over an encounter region. A magnetometer magnetic 
signal value is the magnitude of the projection of a dipole 
magnetic field on the earth magnetic field at the loction of the 
magnetometer. The basis for determining magnetic signal values 
in the program is an expression that involves the magnetic 
dipole moment, the earth magnetic field and the position of the 
magnetometer relative to the dipole. 

The expression can be developed as follows: In the 

rectangular coordinate system that is shown in Figure 1, a 
magnetic dipole is at the origin, the xy-plane is the horizontal 
plane at a representative point in an encounter region, the 
positive y-axis is in the direction of magnetic north, the 
positive x-axis is in the direction of magnetic east and the 
positive z-axis is positive upward. In this rectangular 
coordinate system, the constant earth magnetic field can be 
expressed by Hg = Hg- (j-cos $ + k-sin $) where $ is a 
magnetic dip angle and Hg is a magnetic field magnitude that 
characterizes the earth field in an encounter region. In a 
spherical coordinate system with the origin at the magnetic 
dipole and the polar axis in the direction of the dipole moment, 

= (c • p/r^ ) • ( 2 • r • cos 0 + 0-sin 0) is the magnetic field of the 
dipole at a point with spherical coordinates (r,0,F) . 
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Magnetic 

East 



Figure 1. A unit vector in the direction of the earth field and a 
unit vector in the direction of the dipole field of a 
dipole at the origin are shown at a point that is a 
distance r from the dipole. The unit vector h is 
in the direction of the horizontal component of the 
dipole moment. 
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In this expression p is the magnitude of the dipole moment and 
c is a constant whose value is determined by the choice of 
units. The magnetic signal for a magnetometer that is described 
by the encounter models is Hg = when the magnetometer 

is at the point (r,0,r) . Hg can be expressed in terms of the 
rectangular coordinate system as follows: First, let Tq be a 

unit vector in the direction of the spherical coordinate system 
polar axis. Then the magnetic dipole field 
Hd = (c-p/r^) • (3-r-cos 6 - Tq) 

since r = (i-x + j-y + k-z)/r where r = (x* + y* + z*)5 , 

0 = [ (tq X r) X r]/sin 0 = [(r*ro)-r - (r*r) • rg ]/sin 0, and the 
dot product r«r = cos 0. The unit vector rg can be expressed 
in the rectangular coordinates by noting that rg = p/p and then 
expressing p in rectangular coordinates. To do this, let n 
be the depression angle of p from the xy-plane (the horizontal 
plane) with n positive downward and let a be the direction of 
p relative to magnetic north. Then the magnetic dipole moment 
p = p - (h-cos n - k-sin f2) in terms of the unit vector k and 
the unit vector h = i-sin a + j-cos n which has the direction 
of the horizontal component of p. With these results, the unit 
vector rg = i- (cos n-sin a) + j-(cos n-cos a) - k-sin n and 
Hs = (c- p/r^) • [ (3/r) • (cos 0) • (y • cos (p - z • sin <p) 

“ (cos § • cos n-cos a + sin $-sin n) ] 
where cos 0= (1/r) • (x-cos f2-sin a -f y-cos n-cos a - z-sin n) 
since cos 0 = rTg. As can be seen from this expression, for a 
constant dipole moment magnitude and direction and a constant 
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earth field magnitude and direction, the magnetic signal is only 
a function of the rectangular coordinates of the location of the 
magnetometer relative to the dipole. In the encounter models, 
both of these conditions are satisfied. However, by allowing p, 
n, a, Hg and $ to vary, the expression for Hg is applicable 
to more general encounter models. 
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Appendix 4 . The Anderson Formulation 

In the encounter models, the magnetic signal Hg at a 
sample point in a straight line encounter can be represented in a 
form described by Anderson in Reference 4. For convenience, the 
Anderson formulation is used in the program to determine values 
for Hg. It can be developed as follows: The primed 

rectangular coordinate system that is shown in Figure 2 is 
superimposed on the rectangular coordinate system of Figure 1 so 
that the origin is coincident with the origin of that system. A 
magnetometer is in a straight line encounter with a magnetic 
dipole that is located at the origin of the combined system. The 
combined system moves with the magnetic dipole with the x'-axis 
oriented so that it is parallel to and in the direction of the 
track of the magnetometer relative to the magnetic dipole and the 
z ' -axis oriented so that it is directed toward and passes through 
the CPA on that track. Let 1, m, and n be the direction 
cosines of the dipole moment p and 1]^, mj^ and n^ be the 
direction cosines of the earth magnetic field Hg. Then, the 
unit vector rQ = i' •! + j' *n + k' -m. And, for points on the 
relative track, x' = s', y' =0 and z' = R where R is the 
slant range of the dipole at CPA and where s' is the algebraic 
distance of the magnetometer from CPA on the relative track. (It 
is negative for points before CPA and positive for points after 
CPA.) This implies that magnetometer position vector in the 
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Magnetic 

East 



Figure 2. The dipole is at the origin in the unpriined, primed and 
double primed coordinate systems. In the primed system, the 
coordinates of the magnetometer are (s',0, R) and the 
coordinates of the CPA are (0, 0, R) . 
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moving coordinate system is r = i' • (s'/^) + k' • (R/r) . From 
Appendix 3 , 

Hd = (c-p/r3) • (3 -rcos e - Tq) , 

With r and tq expressed in terms of the primed unit vectors 
and cos © = r*ro = 1- (s'/J^) + it (R/r) > this becomes 

Hd “ (c • p/r^ ) • [ (3/r^) • (1 • s ' + n-R)(i'-s' + k'R) 

- (i' • 1 +j ' • m + k' • n) ] . 

Then, since Hg/Hg = i'-l^ + j'-m]_ + k' • n^ and Hg = Hg*H(j/Hg, 
Hs = (c-p/r^) • [ (2 • 1- li - m-mi - n-ni)-(s')" 

+ 3-(n-li + l-n3^)-s'-R + (2-n-ni - 1-lx - m-mx)-R*]. 
The quantities 

A 2 = 2-1-lx “ m-mx - n-nx 
Ax = 3-(n-lx + l*rix) 

Aq = 2-n-nx “ I'lx “ in-mx, 

are called the Anderson coefficients. With r = [(s')* + R*]^ 

and (3 = s'/R, Hg can now be expressed as follows: 

Hg(P) = (c-p/r 3) -S Aj -Fj (P) 

where the Fj(P) = pl)/(l + P*)^/^ for j =0, 1, 2 are 
called the Anderson functions. This is the Anderson 
formulation . 

To relate the Anderson formulation for Hg to the 
formulation for Hg in Appendix 3, first note that 
1 = ro-i' 
m = ro*j' 
n = ro*k' 

and 
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ll = (He/He) -i' 
mi = (He/He) *j' 
ni = (He/He) -k' . 

Then express Tq, (Hg/Hg) and the unit vectors i' , j' and k' 
in terms of the unit vectors i, j and k and take the 
indicated dot products. From Appendix 3, 

ro = i- (cos n-sin a) + j • (cos Cl-cos a) - k-sin n 

and 

Hq/Hq = j • cos $ + k-sin §. 

To express i' , j' and k' in terms of the unit vectors i, j 
and k, note that the unprimed coordinate system can be 
transformed to the primed coordinate system by two rotations that 
are defined as follows; First, rotate a coordinate system that 
is coincident with the unprimed coordinate system about the 
z-axis through the angle (0 - ir/2) with positive angles clockwise 
(left hand rule) so that its x-axis is parallel to and in the 
direction of the relative track. Next, rotate this system about 
its x-axis through an angle 5 with positive angles 
counterclockwise (right hand rule) so that the positive z-axis 
passes through the CPA. The angle 5 is related to the 
vertical separation Zq of the magnetometer and the dipole and 
the algebraic encounter lateral range L that is positive if the 
dipole is to the left of the relative track. With these sign 
definitions: L = R-sin S and Zq = R-cos 5. After the 

rotation, the auxiliary coordinate system is coincident with the 
primed coordinate system. 
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These transformations can be described in terms of matrix 



equations as follows; Let (x",y",z") be the coordinates of a 
point in the coordinate system that is coincident with the 
auxiliary coordinate system after the first rotation. Then the 
transformation from the unprimed coordinate system to this double 
primed coordinate system is described by the matrix equation 



X" 




sin (/> cos <t> 0 




X 


y" 


= 


-cos 0 sin 0 0 




y 


z" 




0 0 1 




z 



And the transformation from the double primed coordinate system 
to the primed coordinate system is described by the matrix 
equation 



X' 




10 0 




X" 


y 


= 


0 cos 5 sin S 




y" 


z ' 




0 -sin 5 cos 5 




z" 



Taking the product of the rotation matrices in the indicated 
order yields the matrix equation 



X' 




sin 


<t> 






cos 


0 






0 






X 


y’ 


= 


-cos 


5 


COS 


0 


cos 


5 


sin 


0 


sin 


5 




y 


z ' 




sin 


S 


cos 


0 


-sin 


S 


sin 


0 


cos 


5 




z 



which defines the transformation from the unprimed to the primed 
coordinate system. The unprimed vector components of the unit 
vector i' can be found by transforming the coordinates (1,0,0) 
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in the primed system to their corresponding coordinates in the 
unprimed system with the inverse of this matrix and then 
repeating this process for (0,1,0) and (0,0,1) in order to find 
the unprimed unit vector components of j' and k* . However, 
since the inverse transformation matrix is the transpose of this 
matrix, the elements of the row that corresponds to a primed unit 
vector are the magnitudes of the unprimed vectors that are its 
components. Consequently: 

i' = i-sin 0 + j • cos 0 

ji = -i-cos S cos 0 + j • cos S sin 0 + k- sin S 

k' = i-sin S cos 0 - j • sin S sin 0 + k-cos S. 

Then, taking the dot products between r, Tq and these three 

unit vectors as indicate above gives: 

1 = cos n-cos (0 - a) 

m = cos S ■ cos n-sin (0 - a) - sin S • sin 
n = -sin S • cos fi-sin (0 - a) - cos S-sin n 

and 

l] _ = cos $ • cos 0 

m] ^ = cos S • cos $-sin 0 - sin 6 ■ sin $ 

n] ^ = -sin S • cos $-sin 0 - cos 5 • sin 

These are the relations that are used in the program to determine 
values for the Anderson coef f i-cients . 
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Appendix 5. The Encounter Equations of Motion 

In the double primed coordinate system that is defined in 
Appendix 4 , the equations of motion of a magnetometer relative to 
a submarine (dipole) are: 
x” (t) = s * (t) 
y"(t) = -L 
z»(t) = Zq. 

where L is the algebraic encounter lateral range that is 
defined in Appendix 4, Zq is the vertical separation between the 
magnetometer and the submarine and s'(t) is the distance of the 
magnetometer from the CPA on the relative track. With w the 
speed of the magnetometer relative to the submarine and t a 
relative time parameter, s' (t) = w-t . These equations can be 
considered to be the ones used in the program to describe the 
motion of a magnetometer relative to a submarine. There, t is 
determined by t = [j - (m-l)/2]-5t where the index j = 1, 2, 

• • • , m and 5t, a time step, is the time between samples. Note, 

when t = 0, the magnetometer is at the CPA. 

In the coordinate system of Figure 1 in a straight line 
encounter as defined in the encounter models, the equations of 
motion of a magnetometer relative to a submarine can be written 
as follows: 

x(t) = s' (t) -sin 0 + L-cos 0 

y(t) = s' (t) -cos 0 - L-sin 0 

z(t) = Zq, 
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since the transformation from the double primed coordinate system 
to the primed coordinate system is determined by the matrix 
equation 



X 




sin 0 -cos 0 0 




X" 


Y 


= 


cos 0 sin 0 0 




y" 


z 




0 0 1 




z" 



x" = s' (t) and y" = -L. The above equations and the expression 
for Hg in Appendix 3 could have been used in the program to 
evaluate the magnetic signal. In particular, with these 
equations of motion and with the two relations L = R-sin S and 
zq = R-cos 5, the expression for Hg in Appendix 3 can be 
written in terms of n, a, §, <p, S and R so that it is 
identical in appearance to the Anderson formulation for Hg in 
terms of these quantities. The definition of 5 in Appendix 4 
in terms of a counterclockwise rotation results in a definition 
of the algebraic lateral range that is consistent with some that 
have been used elsewhere. 

In the program, the relation w = v - u is the basis for 
determining the relative speed w. In this relation, v is the 
velocity of the magnetometer, u is the velocity of the submarine 
(dipole) and w is the velocity of the magnetometer relative to 
the submarine. This relation implies the following equations: 

Wjj = Vjj - u^/ Wy = Vy - Uy and ~ ~ where the 

coordinates x, y and z refer to a fixed coordinate system 
with the same orientation as that of Figure 1. In the encounters 
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of the models, = v-sin a, Vy = v-cos a and = 0 where o 

is the magnetic course and v is speed of the magnetometer. 

And, Uyr = u sin p, Uy = u-cos P and U 2 = 0 where p is the 
magnetic course and u is the speed of the submarine. The 
relative magnetic course <f> and the relative speed of the 
magnetometer are defined by Wj^ = w-sin 0 and Wy = w-cos 0. In 
the program, 0 and w are determined with a rectangular to 
polar conversion routine where w = ( w^ + Wy) ^ and where 0 is 
determined by sin~^(Wjj/w) and cos~^(Wy/w). 
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Appendix 6. The Submarine Magnetic Dipole 

In the encounter models, the magnitude and direction of a 
submarine's dipole moment p are determined by first determining 
its components in the rectangular coordinate system of Figure 1. 
In that coordinate system, 

Px = (PLP + Pli) ’Sin p + (Ptp + Pti) cos p 
Py = (Plp + Pli) 'COS P - (ptp + Pti) sin p 
Pz = ~(PVP Pvi) 

where P is the submarine's magnetic course, PLp, Ptp and pyp 
are the permanent and Pli^ PtI Pvi induced 

longitudinal, transverse and vertical magnetic dipole moments of 
the submarine. These relations are based on the following sign 
convention: 



PL 


is 


positive 


when 


pl 


is 


directed 


from 


stern to bow. 


pt 


is 


positive 


when 


Pt 


is 


directed 


from 


port to starboard. 


Pv 


is 


positive 


when 


Pv 


is 


directed 


downward • 



The permanent dipole moments are input parameters in the program 
and the induced dipole moments are determined in a way that is 
similar to one that is described in Reference 1. In the 
encounter models, a submarine is a ferromagnetic prolate 
ellipsoid with the major axis the submarine's longitudinal axis 
and the equal minor axes the submarine's transverse and vertical 
axes. And, the induced dipole moments are; 

Pli = ^L'HeL 
Pti ^ ^^T'**eT 
Pvi = ky-Hev 
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where Hql, H^rj. and are the vector components of the earth 

magnetic field in the rectangular coordinate system defined by 
the ellipsoid axes and the submarine magnetic moment sign 
convention. The magnitude of these vector components are: 

HeL Hg-cos # cos (3 
Hg>p = -Hg-cos § sin P 
HgY = Hg-sin 

The earth magnetic field dip angle $ and the earth magnetic 
field magnitude Hg can each be chosen to be input parameters in 
the program or the can be computed by the program as described in 
Appendix 7. (Values for $ and Hg can be found from magnetic 
data charts, for example, see References 5 and References 6.) 

By using the above relations, 

Px ~ Hg-cos $• (kL - krp) - sin p cos P + (pLp-sin p + pfjp-cos P) 

Py = Hg-cos $-(kL-cos*P + kip-sin*P) + (pLp-cos p - p.jip-sin P) 
pz = “(Hg-kysin § + pvp) 

In the encounter models, the values of the permeability 
coefficients kj^, krp and ky are related to submarine 
displacement by the following relations: 
kb = fb’W 
kiji — f ip • W 
kv = fyW 

where W is the submarine displacement in tons and fj^, fj, and 
fy are permeability factors that are determined by a submarine's 
magnetic characteristics. In the program, the units of Hg are 
gamma, the units of p are in oersted-centimeter^, the units of 
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k are oersted-centimeter^/gamma and the units of f are 
oersted-centimeter^/gamma-ton. If the units of Hg were 
gamma, but the units of p were gamma-foot^, then the units of 
k would be foot^ and the units of f would be foot^/ton. To 
convert p in gamma-foot^ to oersted-centimeter^ or k in 
foot^ to oersted-centimeter^/gamma or f in foot^/ton to 
oersted-centimeter^/gamma-ton, divide by 3.53. (The program 
default values for fL, fij> and fy are values from Reference 1 
in foot^/ton that have been divided by 3.53 to give values in 
oersted-centimeter^/gamma-ton. ) To convert p in weber-meter to 
oersted-centimeter^, multiply by [ 1/ (4 • -rr) ] • lO^*^ and to convert 
p in ampere-centimeter^ to oersted-centimeter^ , multiply by 
1- iq2 . 
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Appendix 7. The Earth Magnetic Field and Dip Angle 

An auxiliary magnetic field model is described in this 
appendix. The model is the basis for a default choice for either 
the value of the earth magnetic field magnitude parameter Hg or 
the earth magnetic field dip angle parameter Relative to 

encounter model accuracy, the default values should be addequate 
in most cases. 

In the model, the earth magnetic field is generated by a 
magnetic dipole that is located at the earth's center and the 
earth is a nonmagnetic sphere of radius r^. With Pg the 
magnitude of the dipole moment, the magnitude of the earth field 
at a point is 

Hg = (Hg-Hg)5 = (c-pg/r3) • (3cos*0 + 1)^ 
where © is the polar angle of the point in a spherical 
coordinate system and c is a constant whose value is determined 
by the choice of units. The dipole moment is coincident with and 
in the direction of the polar axis which is directed toward the 
earth's southern hemisphere. In this coordinate system, at any 
point on the surface of the earth: 

He = Hgg- (3 -cos*© + 1) J 

where Hgg is the value of Hg at the magnetic equator which is 
defined by the points on the earth where © = 90°. In terms of 
the dip angle, at any point on the earth's surface; 

Hg = 2Hgo‘ (3-cos*§ + 1)"^. 

This expression can be obtained by noting that f can be 
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defined in terms of the or the r and 0 components of Hg 
the spherical coordinates as follows: 



in 



sin $ = -Hgj-/He = -2 Hgo (cos 0)/Hg 

and 

cos $ = Hg 0 /Hg = Hgo'(sin 0)/Hg. 

Based on these relations, 

[(sin l>)/ 2 ]’ + cosH = (H^q) * 

and 

Hg = HgQ • [ (sin* )/4 + cos*$]“^ = 2 • Hgg • ( 3 • cos * $ + 1 )“^ 
The dip angle is determined from a magnetic latitude for the 
encounter region. 

The magnetic latitude and longitude of the representative 
point of an encounter region can be defined in terms of its 
geographic latitude and longitude by using the following 
transformations: First, convert the latitude and longitude of 

the point to rectangular coordinates in a right-handed coordinate 
system whose origin is at the center of a spherical earth, whose 
positive z-axis passes through its north geographic pole and 
whose positive y-axis passes through the point on its equator 
with latitude 0“ and longitude 0°. Next, rotate a coordinate 
system that is coincident with this system in a clockwise 
direction (left hand rule) about its z-axis so that its positive 
y-axis passes through the point with latitude 0 ° and with 
longitude equal to that of the north magnetic pole. Then, rotate 
the system in a clockwise direction about its x-axis so that its 
positive z-axis passes through the north magnetic pole. Next, 
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transform the rectangular coordinates of the representative point 
in this system to its coordinates in the spherical coordinate 
system that is associated with it. Then, with f2 the polar 
angle of the representative point in this system, the magnetic 
latitude of the representative point is = 90° - n. Since 

tan $ = -2 • cot 0 from above and 0 = 180° - f2, the dip angle # 
is given by the following relation: 

$ = tan“^(2-tan Lj^) . 

In the program, the transformations described above are 
accomplished in part by a rectangular to polar conversion 
routine. In particular, by rotating the final rectangular 
coordinate system about its z-axis so that the x-coordinate is 
zero, the polar angle of the polar coordinates of the 
representative point in the resulting yz-plane determines Ljyj. 

Some values of # and Hq that are listed in Table 1 were 
generated using the program. The latitude and longitude of the 
magnetic pole that are in the program are 76° N and 100° W. 
Values of $ from Reference 6 and of Hg from Reference 7 are 
also listed in Table 1. A comparison of the program values with 
these values gives an indication of the errors inherent in the 
procedure. 
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Encounter 

Latitude 

60°N 
60°N 
60°N 
60°N 
30°N 
30°N 
30 °N 
30°N 
30°N 
30°N 
30°N 
0 ° 

0 ° 

0 ° 

0 ° 

30°S 

30°S 

30°S 

60°S 



Encounter 

Longitude 

180°W 

150°W 

90°W 

30°W 

150°W 

90°W 

60°W 

30°E 

60°E 

90°E 

150°E 

90°W 

30°W 

0 ° 

60°E 

180°W 

90'’W 

90°E 

30°W 



Program 
Dip $ 

74° 

78° 

82° 

75° 

58° 

62° 

59° 

37° 

31° 

30° 

42° 

26° 

9° 
-48 ° 
-25° 
-45° 
-30° 
-62° 
-69° 



Chart #30 
Dip $ 

70° 

73 ° 

84° 

74° 

50° 

60° 

59° 

43 ° 

45° 

44° 

40° 

18° 

-5° 

-25° 

- 20 ° 

-55° 

-34° 

- 66 ° 

-56° 



Program 
Field Hg 

. 63 
. 66 
. 68 
. 64 
, 51 
. 55 
.53 
.41 
. 39 
. 39 
. 43 
. 38 
. 35 
. 35 
. 38 
.44 
. 39 
. 55 
. 60 



Chart #39 
Field Hg 

.53 
. 55 
.61 
.52 
.41 
.51 
.48 
. 42 
.46 
.44 
.40 
. 32 
.28 
.31 
.37 
. 50 
.31 
. 52 
.34 



Table 1. Some corresponding values for the encounter dip angle 
•I and the earth magnetic field intensity Hg in 
gamma. Chart #33 refers to Reference 5 and Chart #39 
refers to Reference 6. 



A small difference between a chart value and a computed value for 
the dip angle $ at some location does not imply a small 
difference between the chart value and the computed value for the 
earth magnetic field intensity at that location. For example, 
note the values for 60 °N latitude and 30 °W longitude. 
(Positive values of <1 indicate the inclination or dip below the 
horizontal of the north seeking end of a dip needle. Negative 
values indicate the inclination below the horizontal of the south 
seeking end.) 
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An extension of the above procedure can be made for finding 
the magnetic variation at a location. However, the values 
generated by using the procedure are generally unsatisfactory. 
Magnetic variation values are charted in Reference 7. 
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Appendix 8. An Alternative Encounter Model 

The alternative model that is referred to in Section III is 
described in more detail in this appendix. In the model, the 
detection range is the slant range R at the CPA for an 
encounter with a specified detection probability (usually .5) is 
defined by: 

R = [c-p/Hs]V3 

where c and p are defined in Appendix 3 and Hg represents 
a minimum detectable average magnetic signal that is defined by; 
Hg = (ORF) -Nm 

where ORF is a signal-to-noise ratio called the operator 
recognition factor and represents the magnetic noise. 

Combining these two relations gives: 

. R = {c-p/[ (ORF) -N^] 

The value for ORF depends both on the specified detection 
probability and on a specified or implied false alarm 
probability. 

The Anderson formulation is consistent with these relations 
in an approximate sense if the average magnetic signal H is 
defined as a root mean square value such that H = (c-p/R^)-k 
and k is an encounter parameter defined by: 
k = (2 [S Aj -Fj (3i) ] 

where the first sum index i = 1, 2, • • • , m and the second sum 

index j = 0, 1, 2. For a particular encounter geometry, k is 
constant and this suggests that the two encounter models could be 
used to determine an average value for k for an encounter 
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region based on average submarine magnetic characteristics. 

Values for both k and H are generated by the program and such 
values can give an indication of the magnitude of the differences 
in detection range estimates that are based on this model and 
either of the other two encounter models . A more detailed 
comparison of these encounter models is described in Reference 8. 
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Appendix 9. An Example of the Program Output 

The program is designed to generate the following 
quantities: encounter parameter values, lateral range function 
values for the crosscorrelation encounter model and for the 
square law encounter model, average magnetic signal values, slant 
range at CPA values, encounter parameter values and magnetic 
signal values. These values can be saved as a program data file 
and/or they can be printed. 

An example of the program's printed output is listed in 
Table 2, Table 3 and Table 4. Figure 3 is a plot of the lateral 
range function values (a lateral range curve) for the square law 
detector that are listed in Table 3. Figure 4 is a plot of the 
magnetic signal values that are listed in Table 4. 
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program file name 
program data file name 
magnetic data file name 
processing data file name 
kinematic data file name 

combined magnetic, processing & kinematic data file name 

encounter latitude (decimal degrees) 

encounter longitude (decimal degrees) 

encounter variation (decimal degrees) 

encounter dip angle (decimal degrees) 

encounter magnetic field intensity (gamma) 

target permanent longitudinal moment (oersted-cm3) 

target permanent transverse moment (oersted-cm3 ) 

target permanent vertical moment (oersted-cm3 ) 

target displacement (tons) 

target longitudinal permeability coefficient 
target transverse permeability coefficient 
target vertical permeability coefficient 
target longitudinal permeability factor 
target transverse permeability factor 
target vertical permeability factor 
sampling period (seconds) 
integration time (seconds) 
adjusted integration time (seconds) 
number of samples per encounter 
magnetometer course (decimal degrees) 
magnetometer speed (knots) 
magnetometer altitude (meters) 
target course (decimal degrees) 
target speed (knots) 
target depth (meters) 

magnetometer relative magnetic course (decimal degrees) 

magnetometer relatiive speed (knots) 

magnetometer-target vertical separation (meters) 

target induced longitudinal dipole moment (oersted-cm3 ) 

target induced transverse dipole moment (oersted-cm3 ) 

target induced vertical dipole moment (oersted-cm3 ) 

target dipole moment (oersted-cm3 ) 

target dipole moment azimuth (decimal degrees) 

target dipole moment depression (decimal degrees) 

distance between samples on the relative track (meters) 

false alarm rate (false alarms per hour) 

false alarm probability 

magnetic noise (gamma) 

maximum lateral range (meters) 

lateral range step (meters) 

number of lateral range function values 



MAD. BAS 
data. mad 
data. mag 
data . pro 
data . kin 
data . mpk 
45 
-45 
-25 

68.41608 

59035.05 
0 

0 

0 

1800 

13140 

2880 

2880 

7.3 

1.6 

1.6 

.5 

20 

20.5 
41 
290 
180 
100 
20 
10 
100 

311.8202 

180.2776 

200 

2.017796E+08 
-4.422566E+07 
1. 58099E+08 
2.601273E+08 
32.63751 
37.42884 
46.37141 
3 

1.666667E-02 
. 35 
500 
20 
51 



Table 2. An example of a encounter parameter values printout. 
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data. mad lateral range function values 



L 

meters 


p(cc) 


p(sl) 


H 

gamma 


R 

meters 


k 


-500 


.2155959 


. 0255824 


. 184247 


538.5165 


2.818488 


-480 


.2721803 


. 029692 


.21034 


520 


2 . 879916 


-460 


. 3464673 


3.582633E-02 


.2410144 


501.5974 


2 . 943744 


-440 


.4415205 


4 . 529992E-02 


.2785504 


483 . 3218 


3 . 009938 


-420 


.5573696 


6. 047795E-02 


. 3231911 


465.1881 


3.078418 


-400 


. 6872897 


. 0856887 


. 3764825 


447 .2136 


3 . 149025 


-380 


.8146549 


. 1287804 


.4404623 


429.4182 


3.221497 


-360 


.9159235 


. 2030898 


.5177304 


411.8252 


3.29544 


-340 


.9745625 


.3272475 


. 6108596 


394 .4617 


3 . 370274 


-320 


.9958122 


. 5141411 


.7233319 


377.3593 


3 .445185 


-300 


.9997211 


.7388175 


.8597726 


360.5551 


3 . 519032 


-280 


.9999951 


. 9177597 


1.025158 


344 . 093 


3 . 590261 


-260 


1 


. 9896971 


1.224529 


328.0244 


3.656775 


-240 


1 


.999714 


1.463644 


312.41 


3.71578 


-220 


1 


.9999992 


1.74981 


297.3214 


3.763606 


-200 


1 


1 


2 . 085394 


282.8427 


3.795517 


-180 


1 


1 


2.470496 


269.0725 


3.805514 


-160 


1 


1 


2.901352 


256.125 


3 . 786234 


-140 


1 


1 


3.359555 


244 . 1311 


3.72898 


-120 


1 


1 


3.819509 


233 . 2381 


3 . 624058 


-100 


1 


1 


4 . 239317 


223 . 6068 


3 .461596 


-80 


1 


1 


4 . 532758 


215.4066 


3.232958 


-60 


1 


1 


4.618943 


208.8061 


2.932828 


-40 


1 


1 


4 .44342 


203.9608 


2.5617 


-20 


1 


1 


3.97739 


200.9975 


2. 128285 


0 


1 


1 


3 .332139 


200 


1 . 651697 


20 


1 


.9999993 


2.540665 


200.9975 


1.165581 


40 


. 9999931 


.908713 


1.76834 


203.9608 


.7394171 


60 


.9933111 


.467368 


1.110188 


208.8061 


. 5637183 


80 


.9998733 


.7874625 


.9039734 


215.4066 


.7782836 


100 


.9999999 


.9856169 


1.086971 


223 - 6068 


1.134662 


120 


1 


.9992399 


1.254863 


233.2381 


1.490748 


140 


1 


.9998478 


1.308583 


244 . 1311 


1.811734 


160 


1 


.9998431 


1.29482 


256.125 


2 . 089946 


180 


1 


.9994995 


1.230946 


269.0725 


2.326072 


200 


1 


.9974561 


1.131703 


282.8427 


2.523877 


220 


1 


. 9872744 


1.017684 


297.3214 


2.688154 


240 


.999999 


.9510013 


.9018685 


312.41 


2.823785 


260 


.9999747 


.8646194 


.7917129 


328.0244 


2.935307 


280 


.9996544 


.7242736 


. 6909791 


344 . 093 


3 . 026728 


300 


.9973592 


.5582789 


. 6011183 


360.5551 


3 . 10149 


320 


.9876337 


.4037092 


. 5222269 


377.3593 


3 . 162495 


340 


. 9608865 


. 2816961 


. 4536729 


394.4617 


3.212155 


360 


.9086407 


. 1948741 


. 3944811 


411.8252 


3 . 252463 


380 


.8301703 


. 1365131 


.3443362 


429.4182 


3 . 285058 


400 


.7333605 


9 . 820431E-02 


.3014853 


447.2136 


3 .311281 


420 


.6299387 


. 0731037 


.2645045 


465.1881 


3.332228 


440 


.5302704 


5. 646636E-02 


. 2325951 


483 .3218 


3.348799 


460 


.4408869 


4 . 522933E-02 


. 205043 


501.5974 


3 .361725 


480 


. 3645446 


3 . 747106E-02 


. 1812233 


520 


3 . 371603 


500 


.3013791 


3 . 199235E-02 


. 1605954 


538.5165 


3 . 378925 



Table 3. An example of a lateral range function values printout. 

The heading for the crosscorrelation values is p(cc) 
and the heading for the square law values is p(sl) . 
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Figure 3. A plot of the square law lateral range function values 
that are listed in Table 3. The horizontal axis is 
encounter horizontal range at CPA in meters. The 
vertical axis is encounter detection probability. 
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data. mad magnetic signal values 
for a lateral range of 40 meters 



relative CPA distance 


magnetic sign< 


s in meters 


Hs in gamma 


-927 .4282 


-1. 619177E-02 


-881.0568 


-1.834315E-02 


-834 . 6854 


-2 . 086144E-02 


-788.314 


-2 . 381637E-02 


-741.9426 


-2 . 728804E-02 


-695.5712 


-.0313641 


-649.1998 


-3. 613118E-02 


-602.8283 


-4 . 165412E-02 


-556.4569 


-4.792992E-02 


-510.0855 


-5.478963E-02 


-463.7141 


-6. 169327E-02 


-417.3427 


-6.730653E-02 


-370.9713 


-.0686387 


-324.5999 


-5.933568E-02 


-278.2285 


-2. 652314E-02 


-231.8571 


.0540608 


-185.4856 


.2215785 


-139.1142 


.5163918 


-92.74282 


.9104893 


-46.37141 


1.184586 


0 


1.003207 


46.37141 


. 3837067 


92.74282 


-.2252113 


139.1142 


-.5323231 


185.4856 


-.5837545 


231.8571 


-.5168925 


278.2285 


-.4202569 


324 . 5999 


-.3304064 


370.9713 


-.2572838 


417.3427 


-.2007385 


463 .7141 


- . 1577988 


510.0855 


-.125285 


556.4569 


-.1005532 


602 .8283 


-8. 158605E-02 


649 . 1998 


-6.689598E-02 


695.5712 


-5. 539922E-02 


741.9426 


-4 . 630719E-02 


788.314 


-3.904362E-02 


834 . 6854 


-3 . 318428E-02 


881.0568 


-.0284142 


927.4282 


-2 . 449729E-02 



Table 4. An example of a magnetic signal values printout. 
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Figure 4. A plot of the magnetic signal values listed in Table 4 
The horizontal axis is relative CPA distance in meters 
The vertical axis is magnetic signal in gamma. 
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Appendix 10. A Program Listing 



10 CIS : CLEAR 

20 EEdUT "Magnetic Anomaly Etetection (MAD) Lateral Range Function Program" 

30 DIMX0(70), R0(70), H(70), HS(70, 150), K(70) , PDCC(70) , PDSL(70) 

40 PI = 3.141592654#: CON = 2 * PI / 360: KON = 1852 / 3600: N$ = "MAD.EAS" 

50 QO = .2316419: Q1 = .31938153#: Q2 = -.356563782#: Q3 = 1.781477937# 

60 Q4 = -1.821255978#: Q5 = 1.330274429# 

70 II = 2.515517: 12 = .802853: 13 = .010328: 14 = 1.432788: 15 = .189269: 16 
= .001308 
80 PRINT 

90 INRJT "generate data or print a program data file (g/p)"; A$ 

100 IF A$ = "G" OR A$ = "g" THEN GOTO 120 

110 IF A$ = "P" OR A$ = "p" THEN GOTO 2960 ELSE GOTO 90 

120 PRINT : A$ = "a" 

130 INEUT "magnetic, processing & kinematic data entry by combined file 
(y/n)"; A$ 

140 IF A$ = "Y" OR A$ = "y" THEN GOTO 2070 

150 IF A$ = "N" OR A$ = "n" THEN GOTO 160 ELSE GOTO 130 

160 INHJT "magnetic data entry by file/keyboard (f/k)"; A$ 

170 IF A$ = "F" OR A$ = "f" THEN GOTO 780 

180 IF A$ = "K" OR A$ = "k" THEN GOTO 190 ELSE GOTO 160 

190 INHJT "latitude in decimal degrees (N +)"; LAT: lATR = LAT * CON 

200 INHJT "longitude in decimal degrees (E +)"; LNG: LNGR = -LNG * CON 

210 INHJT "variation in decimal degrees (E +)"; DEC: DECR = DEC * CON 

220 A$ = "a" 

230 INHJT "input dip angle (y/n)"; A$ 

240 IF A$ = "Y" OR A$ = "y" THEN GOTO 260 

250 IF A$ = "N" OR A$ = "n" THEN GOTO 290 ELSE GOTO 230 

260 INHJT "dip angle in decimal degrees (north magnetic hemisphere +) " ; DIP 
270 DIPR = DIP * OON 
280 GOTO 390 

290 lATPR = 76 * OON: INGPR = 100 * OON 

300 D = SIN(LNGR - INGPR) * COS(IATR) : E = OOS(INGR - INGPR) * COS(IATR) : F = 
SIN(IATR) 

310 X = E: Y = F 
320 GOSUB 4060 

330 T = T - (90 * OON - LATPR) : E = R * SIN(T) : F = R * OOS(T) 

340 X = E: Y = D 
350 GOSUB 4060 
360 X = F: Y = R 
370 GOSUB 4060 

380 DIPR = ATN(2 * TAN(T) ) : DIP = DIPR / OON 
390 A$ = "a" 

400 INPUT "iiput encounter magnetic field intensity (y/n)"; A$ 

410 IF A$ = "Y" OR A$ = "y" THEN GOTO 430 
420 IF A$ = "N" OR A$ = "n" THEN GOTO 450 ELSE GOTO 400 
430 INPUT "encounter magnetic field intensity in gamma"; HE 
440 GOTO 460 

450 HE = 70000! / SQR(3 * COS (DIPR) * COS (DIPR) + 1) 

460 PIM = 0: PIM = 0: PVM = 0 
470 A$ = "a" 
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480 INIur "irpit target peritanent dipole moments (y/n)”; A$ 

490 IF A$ = "Y" OR A$ = "y" THEN GOTO 510 

500 IF A$ = »N" OR A$ = ''n'' THEN GOTO 540 ELSE GOTO 480 

510 INPUT "permanent longitudinal moment in oersted-cm3 (stem-to-bow +)"; PLM 
520 INPUT "permanent transverse moment in oersted-cm3 (port-to-starboard +) " ; 
PIM 

530 INPUT "permanent vertical itonent in oersted-cm3 (downward +) " ; FVM 
540 INPUT "target di^lacement in tons"; WT 
550 A$ = "a" 

560 INPUT "input target permeability coefficients or factors (c/f)"; A$ 

570 IF A$ = "C" OR A$ = "c" THEN GOTO 590 

580 IF A$ = "F" OR A$ = "f" THEN GOTO 630 ELSE GOTO 560 

590 INPUT "longitudinal permeability coefficient in cgs units"; KL 

600 INPUT "transverse permeability coefficient in cgs units"; KT 

610 INPUT "vertical permeability coefficient in cgs xonits"; KV 

620 FL = KL / WTs FT = KT / WT; FV = KV / WT: GOTO 710 

630 INPUT "longitudinal displacement factor in cgs units"; FL 

640 INPUT "transverse displacement factor in cgs units"; FT 

650 INPUT "vertical displacement factor in cgs units"; FV 

660 KL = FL*WT; Kr = FT*Wr: KV = FV*WT 

670 A$ = "a" 

680 INPUT "generate a magnetic data file (y/n)"; A$ 

690 IF A$ = "Y" OR A$ = "y" THEN GOTO 710 

700 IF A$ = "N" OR A$ = "n" THEN GOTO -850 ELSE GOTO 680 

710 INPUT "magnetic data file name"; M$ 

720 OJ ERROR GOTO 730: GOTO 740 

730 RESUME 710 

740 OPEN "O", #1, M$ 

750 WRITE #1, lAT, ING, DEC, DIP, HE, PLM, PIM, PVM, WT, KL, KT, KV, FL, FT, 

FV 

760 CLOSE 
770 GOTO 850 

780 INPUT "magnetic data file name"; M$ 

790 ON ERROR GOTO 800: GOTO 810 

800 RESUME 780 

810 OPEN "I", #1, M$ 

820 INPUT #1, LAT, ING, DEC, DIP, HE, PIM, PIM, PVM, WT, KL, KT, KV, FL, FT, 

FV 

830 CLOSE 

840 DIPR = DIP * CON 
850 A$ = "a" 

860 INPUT "processing data entry by file/keyboard (f/k)"; A$ 

870 IF A$ = "F" OR A$ = "f" THEN GOTO 1450 

880 IF A$ = "K" OR A$ = "k" THEN GOTO 890 EISE GOTO 860 

890 A$ = "a" 

900 INPUT "input sairpling period (y/n)"; A$ 

910 IF A$ = "Y" OR A$ = "y" THEN GOTO 930 

920 IF A$ = "N" OR A$ = "n" THEN GOTO 960 ELSE GOTO 900 

930 INPUT "sanpling period in seconds"; DT 

940 IF DT <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT : GOTO 
930 

950 GOTO 1090 
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960 A$ = "a" 

970 INKfr "ir^ut inaximum magnetic signal frequency (y/n)"; A$ 

980 IF A$ = "Y" OR A$ = "y" THEN GOTO 1000 

990 IF A$ = "N" OR A$ = "n" THEN GOTO 1030 ELSE GOTO 970 

1000 INEUT "maximum magnetic signal frequency in Hertz"; MAXF 

1010 IF MAXF <= 0 THEN PRINT : PE^UT "must be greater than zero": PRINT : GOTO 

1000 

1020 GOTO 1080 

1030 INPUT '"irujiiraum target slant range at CPA in meters"; MINRO 

1040 IF MINRO <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT : 

GOTO 1030 

1050 INPUT "maximum magnetometer relative speed in knots"; MZVXVMK 
1060 MAXVM = MAXVMK * KON: MAXF = 2 * MAXVM / MINRO 

1070 IF MAXVM <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT : 

GOTO 1050 

1080 DT = 1 / (2 * MAXF) : REM low pass filter, Nyquist sanpling rate 
1090 A$ = "a" 

1100 INPUT "irput integration time (y/n)"; A$ 

1110 IF A$ = "Y" OR A$ = "y" THEN GOTO 1130 

1120 IF A$ = "N" OR A$ = "n" THEN GOTO 1210 ELSE GOTO 1100 

1130 INPUT "integration time in seconds"; IT 

1140 IF IT >= ET THEN GOTO 1170 

1150 PRINT : PRINT "IT = " + STR$(IT) + " seconds - minimum = " + STR$(DT) + " 
seconds": PRINT 
1160 GOTO 1130 

1170 NS = 2 * INT(IT / DT / 2) +1: REM adj number of sairples per integration 
time 

1180 IF NS <= 151 THEN GOTO 1340 

1190 PRINT : PRINT "IT = " + STR$(IT) + " seconds - maximum = " + STE^$(150 * 
DT) + " seconds": PRINT 
1200 GOTO 1130 

1210 INPUT ''maximum target slant range at CPA in meters"; MAXRO 

1220 IF MAXRO <= 0 THEN PRINT : PRINT "must be greater than zero": PRINT : 

GOTO 1210 

1230 INPUT ''minimum magnetometer relative speed in knots" ; MINVMK 
1240 MINVM = MINVMK * KOI 

1250 IF MINVM <= 0 THEN PRINT ; PRINT "must be greater than zero": PRINT : 

GOTO 1230 

1260 IT = 2 * MAXRO / MINVM 
1270 IF IT >= DT THEN GOTO 1300 

1280 PRINT : PRINT "IT = " + STR$(IT) + " seconds - minimum = " + STR$(DT) + " 
seconds": PRINT 
1290 GOTO 1100 

1300 NS = 2 * INT(IT / DT / 2) +1: REM adjusted number of samples per 

integration time 

1310 IF NS <= 151 THEN GOTO 1340 

1320 PRINT : PRINT "IT = " + STR$(IT) + " seconds - minimum = " + STR$(150 * 
DT) + " seconds": PRINT 
1330 GOTO 1100 
1340 A$ = "a" 

1350 INPUT "generate a processing data file (y/n)"; A$ 

1360 IF A$ = "Y" OR A$ = "y" THEN GOTO 1380 
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1370 IF A$ = "N” OR A$ = "n" THEN GOTO 1510 EISE GOTO 1350 
1380 INHJT "processing data file name"; P$ 

1390 ON ERROR GOTO 1400: GOTO 1410 
1400 RESUME 1380 
1410 OPEJI "0", #1, P$ 

1420 WRITE #1, DT, IT, NS 

1430 CLOSE 
1440 GOTO 1510 

1450 INR7T "processing data file name"; P$ 

1460 ON ERROR GOTO 1470: GOTO 1480 
1470 RESUME 1450 
1480 OPEN "I", #1, P$ 

1490 INHJT #1, CT, IT, NS 

1500 CLOSE 
1510 A$ = "a" 

1520 INHJT "kinematic data entry by file/keyboard (f/k)"; A$ 

1530 IF A$ = "F" OR A$ = "f" THEN GOTO 1890 

1540 IF A$ = "K" OR A$ = "k" THEN GOTO 1550 EISE GOTO 1510 

1550 INHJT "magnetometer course in decimal degrees (0 if at rest)"; CM 

1560 MCM = (CM - DEC) : lOlR = MCM * OCW: REM magnetometer magnetic course 

1570 INHJT ••magnetometer speed in knots"; VMK 

1580 INHJT ••magnetometer altitude in meters (below 0 is -)"; AM 

1590 INHJT ••target course in decimal degrees^^; CT 

1600 MCT = (CT - DEC) : MCTR = JKT * CON: REM target magnetic course 

1610 INHJT ••target speed in knots"; VTK 

1620 INHJT ••target depth in meters (above 0 is -)••; AT 

1630 IIM = KL * HE * OOS(DIER) * OOS(MCIR) : ITM = -KT * HE * OOS(DIER) * 

SIN (MCTR) 

1640 IVM = KV * HE * SIN(DIRR) 

1650 CMX = (PIM + IIM) * SIN (MCTR) + (PIM + ITM) * COS (MCTR) 

1660 EMY = (PIM + IIM) * OOS(MCTR) - (PIM + ITM) * SIN(MCTR) 

1670 EMV = PVM + IVM 

1680 X = CMX: Y = CMY 

1690 GOSUB 4060 

1700 CMIR = T: X = EMV: Y = R 

1710 GOSUB 4060 

1720 CM = R: CMR = T 

1730 WXK = VMK * SIN(MCMR) - VTK * SIN(MCrR) : WYK = VMK * COS(MCMR) - VTK * 
COS (MCTR) 

1740 Z=AM + AT: REM vertical separation (- for magnetometer below target) 
1750 X = WXK: Y = WYK: GOSUB 4060 

1760 CR = T: WK = R: REM relative magnetometer magnetic course and speed 
1770 CML = CMLR / OON: CM = CMR / CON: C = CR / CON 
1780 A$ = ••a" 

1790 INHJT ••generate a kinematic data file (y/n)"; A$ 

1800 IF A$ = ••Y^* OR A$ = ••y^^ THEN GOTO 1820 

1810 IF A$ = "N^^ OR A$ = •'n^^ THEN GOTO 1950 ELSE GOTO 1790 

1820 INHJT ••kinematic data file name^^; K$ 

1830 ON ERROR GOTO 1840: GOTO 1850 
1840 RESUME 1820 
1850 OPEN "O", #1, K$ 

1860 WRITE #1, CM, VMK, AM, CT, VTK, AT, Z, C, WK, IIM, TIM, IVM, CM, CML, CM 
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1870 CLjOSE 
1880 GOTO 1950 

1890 INHJT "kinematic data file name"; K$ 

1900 ON ERROR GOTO 1910: GOTO 1920 
1910 RESUME 1890 
1920 OPEN "I", #1, K$ 

1930 INPUT #1, CM, VMK, AM, CT, VTK, AT, Z, C, WK, IIM, TIM, IVM, EM, CML, CM 
1940 dOSE 
1950 A$ = "a" 

I960 INPUT "generate a combined magnetic, processing & kinematic data file 
(y/n)"; A$ 

1970 IF A$ = "Y" OR A$ = "y" THEN GOTO 1990 

1980 IF A$ = "N" OR A$ = "n" THEN GOTO 2140 EISE GOTO 1960 

1990 INPUT "combined magnetic, processing & kinematic data file name"; E$ 

2000 ON ERROR GOTO 2010: GOTO 2020 
2010 RESUME 1990 
2020 OPEN "0", #1, E$ 

2030 WRITE #1, lAT, ING, DEC, DIP, HE, PLM, PIM, PVM, WT, KL, KT, KV, FL, FT, 
FV, EXT, IT 

2040 WRITE #1, NS, CM, VMK, AM, CT, VTK, AT, Z, C, WK, IIM, TIM, IVM, EM, OML, 
OM, M$, P$, K$ 

2050 CLOSE 
2060 GCTO 2140 

2070 INPUT "combined magnetic, processing & kinematic data file name"; E$ 

2080 ON ERROR GOTO 2090: GOTO 2100 
2090 RESUME 2070 
2100 OPEN "I", #1, E$ 

2110 INPUT #1, LAT, LNG, DEC, DIP, HE, PIM, PIM, PVM, WT, KL, KT, KV, FL, FT, 
FV, DT, IT 

2120 INPUT #1, NS, CM, VMK, AM, CT, VTK, AT, Z, C, WK, IIM, TIM, IVM, EM, CML, 
CM, M$, P$, K$ 

2130 CLOSE 

2140 CMIR = CML * CON: CMR = CM * CON: DIPR = DIP * CON 
2150 CR = C * CON: W = WK * KDN 

2160 INPUT "required false alarm rate in false alarms per hour"; FAR 
2170 PF = FAR * IT / 3600: REM false alarm probability 
2180 Y = PF: IF PF > .5 THEN Y = 1 - Y: REM inverse normal approximation 
2190 Y = SQR(LOG(l / Y / Y) ) 

2200 Y = Y - (II + Y * (12 + 13 * Y) ) / (1 + Y * (14 + Y * (15 + 16 * Y) ) ) 

2210 IF PF < .5 THEN Y = -Y 
2220 ZP = -Y 

2230 CHI = NS * (1 - 2 / 9 / NS + ZP * SQR(2 / 9 / NS) ) * 3: REM inverse 

chi-square approximation 

2240 INPUT "magnetic noise in gamma"; SIG 

2250 INPUT "maximum lateral range in meters"; IRM 

2260 INPUT "lateral range step in meters"; ST 

2270 IF ST <= IRM THEN GCTO 2290 

2280 PRINT : PRINT "maximum step is " + STR$(IRM) + " meters": PRINT : GOTO 
2260 

2290 NC = 2 * INT(IRM / ST) + 1: REM number of lateral range function values 
2300 IF NC <= 71 THEN GOTO 2320 
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2310 PRINT : PE?INT '*1011111101111 step is " + STR$(I2M / 35) + " meters": PRINT : 
GOTO 2260 

2320 ATT = DT * NS; REM adjijsted integration time 

2330 E6 = W * DT: REM distance between sanples on the relative track 

2340 XO = -(NC - 1) / 2 * ST 

2350 FOR I = 0 TO NC - 1 

2360 X0(I) = XO 

2370 X = XO: Y = Z 

2380 GOSUB 4060 

2390 RO = R; R0(I) = R: RM target slant range at CPA in meters 

2400 DEIR = T: REM target depression angle corrplement at CPA in radians 

2410 IF RO = 0 THEN GOTO 2730: REM zero lateral range and vertical separation 

2420 EMF = CM / 10 / RO * 3: REM dipole moment factor 

2430 L = COS(CMR) * C0S(CR - CMR) 

2440 M = COS (DEIR) * OOS(CMR) * SIN(CR - CmR) - SIN(DEIR) * SIN(CMR) 

2450 N = -SIN(DEIR) * OOS(CMR) * SIN(CR - CMIR) - COS (DEIR) * SIN(C»1R) 

2460 LI = COS(DIPR) * COS(CR) 

2470 Ml = COS (DEIR) * COS(DIPR) * SIN(CR) - SIN(DEIR) * SIN(DIPR) 

2480 N1 = -SIN (DEIR) * OOS(DIIR) * SIN(CR) - COS (DEIR) * SIN(DIER) 

2490 A2 = 2*L*L1-M*M1-N* Nl: REM Anderson Function Coefficient 
2500 Al==3* (N*L1 + L* Nl) : REM Anderson Function Coefficient 
2510 A0 = 2*N*N1-L*L1-M* Ml: REM Anderson Function Coefficient 
2520 SUM = 0: HMAX = 0: HfflN = 0 
2530 FOR J = 0 TO NS - 1 

2540 S = (J - (NS - 1) / 2) * IDS: BA = S / RO: REM Anderson Function Argument 

2550 AF = 1/ (1 + BA* BA) * 2.5: REM Anderson Function Factor 

2560 HSF = (A2 *BA*BA + A1*BA + AO) * AF: REM magnetic signal factor 

2570 HS(I, J) = CMF * HSF: REM magnetic signal value 

2580 IF HS(I, J) > HMAX THEN HMAX = HS(I, J) 

2590 IF HS(I, J) < HMIN THEN HMIN = HS(I, J) 

2600 SUM = SUM + HSF * HSF 
2610 NEXT J 

2620 H(I) = HMAX - HMIN 
2630 K(I) = SQR(SUM) 

2640 W = -ZP + CMF * SQR(SUM) / SIG 

2650 LAM = CMF * CMF * SUM / (SIG * SIG) : AA = NS + lAM: BB = 1 + LAM / (NS + 
LAM) 

2660 ZN = -SQR(2 * CHI / BB) + SQR(2 *AA/BB-1): X1 = W 

2670 GOSUB 4120 

2680 IF Y1 > 1 THEN Y1 = 1 

2690 PDCC(I) = Yl: XI = ZN 

2700 GOSUB 4120 

2710 IF Yl > 1 THEN Yl = 1 

2720 PDSL (I) = Yl 

2730 XO = XO + ST 

2740 NEXT I 

2750 A$ = "a" 

2760 INPUT "generate a program data file (y/n) " ; A$ 

2770 IF A$ = "Y" OR A$ = "y" THEN GOTO 2790 

2780 IF A$ = "N" OR A$ = "n" THEN GOTO 3120 ELSE GOTO 2760 

2790 INPUT "program data file name"; D$ 

2800 ON ERROR GOTO 2810: GOTO 2820 
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2810 RESUME 2790 
2820 OPEN "0”, #1, D$ 

2830 WRITE #1, LAT, ING, DEC, DIP, HE, PIU, PIM, PVM, WT, KL, KT, KV, 
F7, DT, IT, ATT 

2840 WRITE #1, NS, CM, VM, AM, CT, VT, AT, IIM, TIM, IVM, Z, C, W, EM 
(M, EAR, PF, SIG, ST 

2850 WRITE #1, IRM, DS, NC, M$, P$, K$, E$ 

2860 FOR I = 0 TO NC - 1 

2870 WRITE #1, XO(I) , PDCC(I) , PDSL(I) , K(I) , H(I) , RO(I) 

2880 NEXT I 

2890 FOR I = 0 TO NC - 1 
2900 FOR J = 0 TO NS - 1 
2910 WRITE #1, HS(I, J) 

2920 NEXT J 
2930 NEXT I 
2940 CIDSE 
2950 GOTO 3120 

2960 INPUT "program data file name"; D$ 

2970 ON ERROR GOTO 2980: GOTO 2990 
2980 RESUME 2960 
2990 OPEN "I", #1, D$ 

3000 INPUT #1, lAT, ING, DEC, DIP, HE, PIN, PIM, PVM, WT, KL, KT, KV, 
FV, nr, IT, AIT 

3010 INPUT #1, NS, CM, VM, AM, CT, VT, AT, IIM, TIM, IVM, Z, C, W, EM 
CM, FAR, PF, SIG, ST 

3020 INPUT #1, IRM, DS, NC, M$, P$, K$, E$ 

3030 FOR I = 0 TO NC - 1 

3040 INPUT #1, X0(I), PDCC(I), PDSL(I) , K(I) , H(I) , R0(I) 

3050 NEXT I 

3060 FOR I = 0 TO NC - 1 

3070 FOR J = 0 TO NS - 1 

3080 INPUT #1, HS(I, J) 

3090 NEXT J 
3100 NEXT I 
3110 CLOSE 

3120 PRINT : A$ = "a" 

3130 INPUT "print encounter parameter values (y/n)"; A$ 

3140 IF A$ = "Y" OR A$ = "y" THEN GOTO 3160 

3150 IF A$ = "N" OR A$ = "n" THEN GOTO 3670 EtSE GOTO 3130 

3160 LPRINT 

3170 LPRINT "program file name 
3180 LPRINT "program data file name 
3190 LPRINT "magnetic data file name 
3200 LPRINT "processing data file name 
3210 LPRINT "kinematic data file name 

3220 LPRINT "combined magnetic, processing & kinematic data file name 
3230 LPRINT "encounter latitude (decimal degrees) ": 

lAT 

3240 LPRINT "encounter longitude (decimal degrees) ", 

ING 

3250 LPRINT "encounter variation (decimal degrees) ", 

DEC 



FL, FT, 
CML, 



FL, FT, 
CML, 



" + N$ 

" + D$ 

" + M$ 

" + P$ 

" + K$ 

" + E$ 
SPC(2) ; 

SPC(2) ; 

SPC(2) ; 
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SPC(2) 



3260 LERINT 
DIP 

3270 LPRINT 
HE 

3280 LPRINT 
PIM 

3290 LPRINT 
PIM 

3300 LPRINT 
PVM 

3310 LERINT 
WT 

3320 LPRINT 
KL 

3330 LPRINT 
KT 

3340 LPRINT 
KV 

3350 LERINT 
FL 

3360 LPRINT 
FT 

3370 LERINT 
FV 

3380 LPFONT 
DT 

3390 LPRINT 

rr 

3400 LERINT 
ATT 

3410 LPRINT 
NS 

3420 LPRINT 
CM 

3430 LERINT 
VMK 

3440 LPRINT 
AM 

3450 LPRINT 
CT 

3460 LPRINT 
VTK 

3470 LPRINT 
AT 

3480 LPRINT 
C 

3490 LERINT 
WK 

3500 LPRINT 
Z 

3510 LPRINT 
IIM 



"encounter dip angle (decimal degrees) " 
"encounter magnetic field intensity (gamma) " 
"target permanent longitudinal moment (oersted-cm3) " 
"target permanent transverse moment (oersted-cm3) " 
"target permanent vertical moment (oersted-cm3) " 
"target displacement (tons) " 
"target longitudinal permeability coefficient " 
"target transverse permeabilil^ coefficient " 
"target vertical permeabilily coefficient " 
"target longitudinal permeability factor " 
"target transverse permeabilil^ factor " 
"target vertical permeability factor " 
"sampling period (seconds) " 
"integration time (seconds) " 
"adjusted integration time (seconds) " 
"number of samples per encounter " 
"magnetometer course (decimal degrees) " 
"magnetometer speed (knots) " 
'•magnetometer altitude (meters) " 
••target course (decimal degrees) " 
'•target speed (knots) '• 
"target depth (meters) " 
••magnetometer relative magnetic course (decimal degrees) '• 
••magnetometer relative speed (knots) ••; 
••magnetometer-target vertical separation (meters) " 
"target induced longitudinal dipole moment (oersted-cm3) " 



SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2 ) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 

SPC(2) ; 
; SPC(2); 
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3520 LPRINT "tcLcget induced transverse dipole moment (oersted-cm3) " 

riM 

3530 LEKENT "target induced vertical dipole moment (oersted-cm3) " 
IVM 

3540 LFRINT "tctrget dipole moment (oersted-cm3) " 

CM 

3550 LPRINT "target dipole moment azimuth (decimal degrees) " 

CML 

3560 LPRINT "target dipole moment depression (decimal degrees) " 

CM 

3570 LPRINT "distance between samples on the relative track (meters) " 
DS 

3580 LPRINT "false alarm rate (false alarms per hour) " 

FAR 

3590 LPRINT "false alarm probability " 

PF 

3600 LFRINT "magnetic noise (gamma) " 

SIG 

3610 LFRINT "maximum lateral range (meters) " 

IRM 

3620 LFRINT "lateral range step (meters) " 

ST 

3630 LPRINT "number of lateral range function values " 

NC 

3640 FOR I = 0 TO 15 
3650 LFRINT 
3660 NEXT I 
3670 PRINT ; A$ = "a" 

3680 INPUT "print lateral range function values (y/n) " ; A$ 

3690 IF A$ = "Y" OR A$ = "y" THEN GOTO 3710 

3700 IF A$ = "N" OR A$ = "n" THEN GOTO 3790 ELSE GOTO 3680 

3710 LFRINT E)$; " lateral range function values" 

3720 LFRINT : LPRINT 

3730 LPRINT "L p(CC) p(sl) H R 

k" 



SPC(2) ? 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPG (2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ; 
SPC(2) ? 



3740 LFRINT "meters gamma meters" 

3750 LFRINT 

3760 FOR I = 0 TO NC - 1 

3770 LFRINT X0(I) ; TAB(IO) ; PDCC(I) ; TAB(24) ; PDSL(I) ; TAB(38) ? H(I) ; TAB(52) ; 
R0(I); TAB(70); K(I) 

3780 NEXT I 

3790 PRINT : A$ = "a" 

3800 INFUT "print magnetic signal values (y/n)"; A$ 

3810 IF A$ = "Y" OR A$ = "y" THEN GOTO 3830 

3820 IF A$ = "N" OR A$ = "n" THEN GOTO 4010 ELSE GOTO 3800 

3830 PRINT : PRINT "identify signal by encounter lateral range index" 

3840 PRINT "lateral range equals the index times " + STR$(ST) + " meters" 

3850 PRINT "index values; + STR$((NC - 1) / 2) + " to " + STR$((NC - 1) / 

2) 

3860 INFUT "lateral range index"; K: I = K + (NC - 1) / 2 

3870 LPRINT : LFRINT : LFRINT 

3880 LPRINT E)$; " magnetic signal values" 
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3890 LEEUNT "for a lateral range of K * ST; " meters" 

3900 LPRINT 

3910 LERINT "relative CPA distance"; TAB(35) ; "magnetic signal" 

3920 LFRINT "s in meters"; TAB(35) ; "Hs in gamma" 

3930 LPRINT 

3940 FOR J = 0 TO NS - 1 

3950 IPRINT (J - (NS - 1) / 2) * DS; TAB(35) ; HS(I, J) 

3960 NEXT J 

3970 PRINT : A$ = "a" 

3980 INPUT "print magnetic signal values for a different lateral range (y/n)" 
A$ 

3990 IF A$ = "Y" OR A$ = "y" THEN GOTO 3830 

4000 IF A$ = "N" OR A$ = "n" THEN GOTO 4010 EISE GOTO 3980 

4010 PRINT ; A$ = "a" 

4020 INPUT "continue to use the program (y/n)"; A$ 

4030 IF A$ = "Y" OR A$ = "y" THEN GOTO 10 

4040 IF A$ = "N" OR A$ = "n" THEN GOTO 4050 ELSE GOTO 4020 

4050 END 

4060 R = SQR(X * X + Y * Y) : REM rectangular to polar conversion 
4070 IF R = 0 THEN T = 0: RETURN 

4080 IF ABS(X / R) = 1 THEN Q = SGN(X) * (PI / 2) ELSE Q = AIN(X / R / SQR(1 
X * X / R / R) ) 

4090 IF ABS(Y / R) = 1 THEN T = (PI / 2) * (1 - SQT(Y) ) ELSE T = (PI / 2) - 
AIN(Y / R / SQR(1 - Y*Y/R/R)) 

4100 IF Q < 0 THEN T = 2 * PI - T 
4110 REIURN 

4120 Y1 = XI: IF XI < 0 THEN Y1 = -Yl; REM normal approximation 
4130 G = 1 / (1 + QO * Yl) 

4140 Yl = EXP(-Y1 * Yl / 2) / SQR(2 * PI) * G * (Q1 + G * (Q2 + G * (Q3 + G * 
(Q4 + G * Q5) ) ) ) 

4150 IF XI >= 0 THEN Yl = 1 - Yl 
4160 RETURN 
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